############################################################################
### Script for Centers of Gravity: Regional Powers, Democracy, and Trade ###
############################################################################

## Note: code for intial directory setup and models is included here. Coding of variables, as well as creation of table and figures, ##
## is presented in additional script files. ##

library(countrycode)
library(car)
library(foreign)
library(plyr)
library(plm)
library(DataCombine)
library(lme4)
library(effects)
library(ggplot2)
library(stargazer)

## Set working directory for files ##
setwd("/Users/timothypeterson/Documents/Dropbox/Regional Hegemony/II acceptance/II replication")

# Create directories for output
dir.create(path=paste0(getwd(), "/figures"))
dir.create(path=paste0(getwd(), "/tables"))



##################################
### Run code for country-years ###
##################################

source("2 country-year code.R")



###################
### Main models ###
###################

# cyear intra-regin trade - excl. hegemon
model1 <- plm(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
               data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
               na.action = na.omit, method="nr")
summary(model1)


## redo M1 with manual fe (for figure 1) ##
model1.man <- lm(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum + fyear + fccode1,
               data = final.data[final.data$hegemon == 0,])
summary(model1.man)


# cyear intra-region incl hegemon
model2 <- plm(log.intratrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model2)


# cyear inter-region excl hegemon
model3 <- plm(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
               data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
               na.action = na.omit, method="nr")
summary(model3)

## redo model 3 with manual fe (for figure 2) ##
model3.man <- lm(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum + fyear + fccode1,
               data = final.data[final.data$hegemon == 0,])
summary(model3.man)


# cyear inter-region incl hegemon
model4 <- plm(log.intertrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model4)


# Intra-region proportion - nh
model5 <- plm(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
               data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model5)


## Redo model 5 with manual fe (for figure 3) ##
model5.man <- lm(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum + fyear + fccode1,
               data = final.data[final.data$hegemon == 0,])
summary(model5.man)


# Intra-region proportion - all
model6 <- plm(pertrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model6)

## Calculations for text ##
## logged trade values necessitate more complex calculation ##

100*(exp(0.31+(10*.25)) - 1)

100*(exp(0.31+(-10*.25)) - 1)

100*(exp(0.31+(0*.25)) - 1)

100*(exp(-1.17+(10*.03)) - 1)

100*(exp(-1.17+(-10*.03)) - 1)

# Regional trade proportion is non-logged, so the only calculation is for the interaction
.19 + (-10*0.01)

.19 + (10*0.01)

######################
### Appendix stuff ###
######################

## random effects at region and country level ##

# RE version of m1
model1.re <- lmer(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum + count + fyear
                   + (1 | lemkeregion) + (1 | ccode1), data = final.data[final.data$hegemon == 0,]) 
summary(model1.re)

# RE version of m2
model2.re <- lmer(log.intratrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum + count + fyear
                  + (1 | lemkeregion) + (1 | ccode1), data = final.data[final.data$hegemon == 0,]) 
summary(model2.re)

# RE version of m3
model3.re <- lmer(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum +  count + fyear
                  +  (1 | lemkeregion) + (1 | ccode1), data = final.data[final.data$hegemon == 0,]) 
summary(model3.re)

# RE version of m4
model4.re <- lmer(log.intertrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum +  count + fyear
                  +  (1 | lemkeregion) + (1 | ccode1), data = final.data[final.data$hegemon == 0,]) 
summary(model4.re)

# RE version of m5
model5.re <- lmer(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum +  count + fyear
                 + (1 | lemkeregion) + (1 | ccode1), data = final.data[final.data$hegemon == 0,])
summary(model5.re)

# RE version of m6
model6.re <- lmer(pertrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum +  count + fyear
                  + (1 | lemkeregion) + (1 | ccode1), data = final.data[final.data$hegemon == 0,])
summary(model6.re)


## Check results of main models by region ##

# Model 1; Remove Western Europe
model1.n4 <- plm(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
                 data = final.data[final.data$hegemon == 0 & final.data$lemkeregion != 4,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                 na.action = na.omit, method="nr")
summary(model1.n4)

# Model 1;Remove Western Europe and North America
model1.n4.1 <- plm(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
                 data = final.data[final.data$hegemon == 0 & final.data$lemkeregion != 4 & final.data$lemkeregion != 1,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                 na.action = na.omit, method="nr")
summary(model1.n4.1)


# Model 3; Remove Western Europe
model3.n4 <- plm(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0 & final.data$lemkeregion != 4,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model3.n4)

# Model 3;Remove Western Europe and North America
model3.n4.1 <- plm(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
                 data = final.data[final.data$hegemon == 0 & final.data$lemkeregion != 4 & final.data$lemkeregion != 1,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                 na.action = na.omit, method="nr")
summary(model3.n4.1)

# Model 5; Remove Western Europe
model5.n4 <- plm(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0 & final.data$lemkeregion != 4,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model5.n4)

# Model 5;Remove Western Europe and North America
model5.n4.1 <- plm(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
                 data = final.data[final.data$hegemon == 0 & final.data$lemkeregion != 4 & final.data$lemkeregion != 1,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                 na.action = na.omit, method="nr")
summary(model5.n4.1)




## Check results of main models by time period ##

# Model 1; Cold War
model1.60 <- plm(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0 & final.data$year > 1959 & final.data$year < 1989,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model1.60)

# Model 1; post-Cold War
model1.80 <- plm(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
                 data = final.data[final.data$hegemon == 0 & final.data$year > 1988,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                 na.action = na.omit, method="nr")
summary(model1.80)

# Model 3; Cold War
model3.60 <- plm(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0 & final.data$year > 1959 & final.data$year < 1989,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model3.60)

# Model 3; post-Cold War
model3.80 <- plm(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
                 data = final.data[final.data$hegemon == 0 & final.data$year > 1988,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                 na.action = na.omit, method="nr")
summary(model3.80)

# Model 5; Cold War
model5.60 <- plm(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
              data = final.data[final.data$hegemon == 0 & final.data$year > 1959 & final.data$year < 1989,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
              na.action = na.omit, method="nr")
summary(model5.60)

# Model 5; post-Cold War
model5.80 <- plm(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2 + conflict.dum,
                 data = final.data[final.data$hegemon == 0 & final.data$year > 1988,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                 na.action = na.omit, method="nr")
summary(model5.80)

## Models excluding conflict dummy

# cyear intra-regin trade - excl. hegemon
model1nc <- plm(log.intratrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2,
                data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                na.action = na.omit, method="nr")
summary(model1nc)

# cyear intra-region incl hegemon
model2nc <- plm(log.intratrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2,
                data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                na.action = na.omit, method="nr")
summary(model2nc)


# cyear inter-region excl hegemon
model3nc <- plm(log.intertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2,
                data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                na.action = na.omit, method="nr")
summary(model3nc)

# cyear inter-region incl hegemon
model4nc <- plm(log.intertrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2,
                data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                na.action = na.omit, method="nr")
summary(model4nc)

# Intra-region proportion - nh
model5nc <- plm(pertrade.nh1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2,
                data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                na.action = na.omit, method="nr")
summary(model5nc)

# Intra-region proportion - all
model6nc <- plm(pertrade.all1 ~ NHHI * hegpolity + totcinc + log.gdp + log.pop + p_polity2,
                data = final.data[final.data$hegemon == 0,], model = "within", effect = "twoways", index = c("ccode1", "year"), 
                na.action = na.omit, method="nr")
summary(model6nc)



###############################
### Run code for dyad-years ###
###############################

source("3 dyad-year code.R")


## Intra-region trdae dyadic models ##

# Intra-region dyadic model, region country, and dyad random effects
model1dr <- lmer(log.exp1to21 ~ NHHI1 * hegpolity1 + totcinc1 + count1 + log.gdp1 + log.gdp2
                + log.distance + log.arabrat + log.poprat + log.caprat +
                  (1 | lemkeregion1) + (1 | fdyad) + (1 | ccode1) + (1 | ccode2), data = sameregion.nh ) 
summary(model1dr)

# Intra-region dyadic model, region, dyad, and country-year effects - only works with intra-region subset 
# (R crashes from size of inter-region version)
model1dry <- lmer(log.exp1to21 ~ NHHI1 * hegpolity1 + totcinc1 + count1 + log.gdp1 + log.gdp2
                 + log.distance + log.arabrat + log.poprat + log.caprat +
                   (1 | lemkeregion1) + (1 | fdyad) + (1 | ccode1:fyear) + (1 | ccode2:fyear), data = sameregion.nh ) 
summary(model1dry)

# Manual two-way fe - only works with intra-region subset (R crashes from size of inter-region version)
model1dfe <- lm(log.exp1to21 ~ NHHI1 * hegpolity1 + totcinc1 + log.gdp1 + log.gdp2
            + log.arabrat + log.poprat + log.caprat + fyear + fdyad, data = sameregion.nh)
summary(model1dfe)


# Inter-region dyadic model, region country, and dyad random effects
model2dr <- lmer(log.exp1to21 ~ NHHI1 * hegpolity1 + totcinc1 + count1 + NHHI2 * hegpolity2 + totcinc2 + count2
              + log.gdp1 + log.gdp2 + log.distance + log.arabrat + log.poprat + log.caprat +
                (1 | lemkeregion1) +  (1 | lemkeregion2) + (1 | fdyad) + (1 | ccode1) + (1 | ccode2), data = difregion.nh)
summary(model2dr)


## Quasipoisson models ##

# QP version of intra-region dyadic model
modelP1 <- glm(log.exp1to21 ~ NHHI1 * hegpolity1 + totcinc1 + count1 + log.gdp1 + log.gdp2
               + log.distance + log.arabrat + log.poprat + log.caprat, data = sameregion.nh,
               family = quasipoisson(link="log"))
summary(modelP1)

# QP version of inter-region dyadic model
modelP2 <- glm(log.exp1to21 ~ NHHI1 * hegpolity1 + totcinc1 + count1 + NHHI2 * hegpolity2 + totcinc2 + count2 +log.gdp1 + log.gdp2
                + log.distance + log.arabrat + log.poprat + log.caprat, data = difregion.nh,
                family = quasipoisson(link="log"))
summary(modelP2)




#########################
### Table and Figures ###
#########################

# This code produces all tables and figures used in the article and appendix.
source("3 tables and figures.R")

